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In this paper we describe the design and implementation of TARANG, a pseudospectral code 
to simulate turbulent flows in fluids, magnetohydrodynamics (MHD), convection, passive scalar, 
etc. We use the object-oriented features of C++ to abstract operations involved in the simulation. 
TARANG has been validated and used for solving problems in convection and MHD. 



INTRODUCTION 

Turbulence is one of the most challenging and unsolved 
problems of physics. Theoretical or exact analysis of 
turbulence has been rare due to the complex nonlineari- 
ties present in the equations. Therefore, major attempts 
to understand turbulence has been through experiments 
and numerical simulations. Over time very powerful 
supercomputers and software tools have emerged, that 
have propelled the computational capabilities of turbu- 
lent flows to unimaginable heights. 

Some of the popular schemes to solve fluid flows are 
finite difference, finite element, finite volume, spectral el- 
ements, pseudo-spectral, vortex method, etc. The pseudo- 
spectral method pQ is most accurate among them, and 
it employed for studying small-scale turbulence. In the 
present paper we will describe the design and capability 
of a pseudo-spectral code named "TARANG" , which has 
been developed by our group. TARANG is a Sanskrit 
word that means "waves" . 

TARANG is an object-oriented parallel pseudo- 
spectral code that can simulate flows in fluids, mag- 
nctofluids, convection, magnetoconvection etc. The con- 
vective flow module can also be used to solve Rayleigh- 
Taylor instability and turbulence, stratified flows, non- 
Boussinesq convection, and passive scalars, etc. We de- 
signed TARANG in an object-oriented fashion for gener- 
ality and easy adoptability to solve varied problems. We 
make use of fast libraries, FFTW (Fastest Fourier Trans- 
form in the West) and blitz++ for efficiency. The code is 
neatly segregated into libraries and src directory. The 
fluid, magnetohydrodymics (MHD), convection solvers 
etc. make use of the libraries, and they are arranged 
in the src directory. 

We will describe the design issues and various features 
of TARANG in the following sections. 



DESIGN AND IMPLEMENTATION ISSUES OF 
TARANG 

In TARANG we use the object-oriented features of 
C++ to design general purpose libraries to create solvers 
for the incompressible fluid flows. For example, a library 
function Computejnlin computes (u • V)u for all basis 



functions. A solver containing many complex features 
and boundary conditions, and more field variables are 
easily constructed using these functions. Main features 
including the class structures of TARANG are described 
below. 



External libraries 

The handling of arrays and mathematical functions is 
computationally slow in C++. Fortunately, several effi- 
cient C++ libraries are available to perform the above 
tasks. We chose blitz++ for TARANG since it handles 
multidimensional arrays in a nice and succinct manner. 
The other libraries with similar functions are boost, ndar- 
rays, and eigen, but presently we continue our develop- 
ment with blitzH — K 

Pseudospectral codes use FFT (Fast Fourier Trans- 
forms) heavily. In a typical code, approximately 80% 
of the computational time is spent of FFT. FFTW is the 
most popular and the most efficient parallel FFT library 
available today, therefore we use FFTW in our code. 

Basis functions 

A pseudo-spectral code uses basis functions to ex- 
pand the real-space functions. The choice of the basis 
functions depend critically on the boundary conditions, 
exp (ik • x) , where k, x are the wavenumber and real- 
space coordinates respectively, is the natural choice for 
periodic boundary conditions. Convection, channel flows 
etc. however involves walls. All the components of the ve- 
locity fields at the wall must vanish for the no-slip bound- 
ary condition. Chebyshev and Legandre polynomials are 
used to expand such functions. For the free-slip bound- 
ary condition, the velocity field perpendicular to the wall, 
and the perpendicular gradient of the horizontal velocity 
components are zero. Sine and cosine functions are ob- 
vious and simple choices for such simulations. Spherical 
harmonics are natural choice for spherical simulations. 

TARANG focusses on basis-independent libraries, so 
we have designed the basis functions in a modular 
fashion. At present TARANG has FOUR and SCFT 
(sin/cos-Fourier) basis functions that can simulate flows 
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in a box geometry under periodic boundary conditions 
(all directions) and free-slip boundary conditions along 
x, and periodic along y and z directions. The SCFT 
basis functions with appropriate modifications has been 
also used to simulate flows with free-slip boundary con- 
ditions along all the directions. The no-slip boundary 
conditions has been successfully tested for channel flow, 
but this module will be integrated with TARANG soon. 
The spherical geometry and cylindrical geometry will be 
implemented in future. 

Primary functions related to the basis functions are 
forward-transform (transformation from the real space to 
Fourier space), inverse-transform (transformation from 
the Fourier space to the real space), computation of en- 
ergy spectrum etc. We use FFTW for majority of trans- 
form operations. 

Objects of TARANG 

The class structure of TARANG is illustrated in Fig. [I] 
The class IncFluid contains the Incompressible velocity 
field and the associate solvers. This is the final class that 
inherits more classes. We describe the main features of 
the classes below. 



Class CVF: Complex Vector Field 

The class CVF stands for Complex Vector Field. It 
contains three dynamic arrays associated with the three 
components of the velocity or magnetic fields in the 
Fourier space. As mandated by FFTW, the size of 
each array is N%, N 2 , N3/2 + 1 that spans wavenumbers 
{ki,k 2 ,k 3 ) = [-JVi/2 : N 1 /2,-N 2 /2 : N 2 /2,0 : N 3 /2] in 
FOUR basis, and (k u k 2 ,k 3 ) = [0 : N 1: -N 2 /2 : N 2 /2, : 
N3/2] in SCFT basis. These arrays contains complex 
numbers that represent the Fourier amplitude of the vec- 
tor field. The arrays are created dynamically at the run- 





CSF 












CVF RVF 


NLIN 


EnergyTr 








IncVF 


Time 













IncFluid 



time. 

Forward and inverse transforms, and input/output of 
the vector fields are some of the main functions of the 
class CVF. 



Class RVF: Real Vector Field 

The class RVF stands for Real Vector Field, and it con- 
tains three dynamic arrays to represent the vector fields 
in the real space. We still create complex arrays of the 
size Ni,N 2 , A3/2 + 1 for ease of FFTW and blitz++ op- 
erations; here the real and imaginary parts of a complex 
number represent two adjacent points in the real space. 



Class CSF: Complex Scalar Field 

The class CSF stands for Complex Scalar Field. It 
contains a dynamic arrays associated with a scalar field, 
e.g., temperature in the Fourier space. The indexing of 
the array is similar to that of CVF. 



Class RSF: Real Scalar Field 

The class RSF stands for Real Scalar Field, and it con- 
tains a dynamic array associate with a scalar. The array 
features are same as that for RVF. 

The above four classes reside in directory named fields. 

Class IncVF: Incompressible Vector Field 

The class IncVF, acronym for Incompressible Vector 
Field, contains most crucial functions of the solver. In- 
cVF inherits classes CVF, RVF, NLIN, EnergyTr. The 
classes CVF, RVF contain the velocity field in the Fourier 
and real space respectively. The class NLIN contains 
three arrays for storing the nonlinear term djUjUi, where 
the symbol ^represents the Fourier transform. In addi- 
tion, NLIN inherits CSF, whose array is used for storing 
the pressure field. The class EnergyTr contains function 
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FIG. 1. The class structure of IncFluid (Incompressible Flu- 
ids), which is the final class of TARANG. 



FIG. 2. The class structure of IncSF (Incompressible scalar 
field). 
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for computing energy flux, shell-to-shell energy transfer 
etc. 

The class IncVF also contains three arrays Force-i 
to store the force fields, and two array *VF_temp, 
*VF_temp2 to save temporary fields. It also has an array 
*VF_temp_r that is used for storing temporary arrays in 
real space. 

Class IncSF: Incompressible Scalar Field 

The class IncSF is used for a scalar field accompanying 
incompressible velocity field. For example, in Rayleigh 
Benard convection this class is used to represent the tem- 
perature field. IncSF inherits a CSF and a RSF to store 
the scalar field in the Fourier and real space respectively 
(see Fig. [2J. It also contains arrays nlin, Force, and 
SFjtemp to store nonlinear term u • VT, forcing, and 
temporary array. 

Corapute_nlin() : a class function 

The class IncVF has many functions. However, Com- 
putc-nlin is one of the most important functions of this 
class. We will describe this function as an illustration of 
TARANG function: 

void IncVF :: Compute_nlin () 
{ 

*Vlr = *V1; 
*V2r = *V2; 
*V3r = *V3; 

// Inverse transform of *Vir using 

// *VF_temp_r as temporary array 

RV_ Inver se_t r ansf orm ( * VF_t emp_r ) ; 

// Vr[i] -> Vr[i]~2 stored in nlin[i] 

Compute_RSprod_diag () ; 

// nlin[i]= Di T [Vr [i] ~2] ; 

// T = Forward transform 

// Di = deriv ativ e along i-th dim 

NLIN_diag_Forward_transf orm_der ivat ive: 

(* VF_temp_r ) ; 
// Vr[i] = Vr [i] * Vr [j] 
Compute_RSprod_of f diag () ; 
// Vr[i] = T (Vr [i] *Vr [j] ) 
RV_For war d_tr ansf orm_RSprod (* 

VF_temp_r ) ; 
// nlinlil = Dj[T(Uj * Ui ) ] 
Derivative_RSprod_VV () ; 

} 

Listing 1. Compute_nlin 

The comments above the C++ statements explain the 
logic of the functions. Related functions compute the 
nonlinear term in the presence of scalar field and another 
vector field. 



• void Compute Jilin_scalar(IncSF& T): nlin_i = 
djUjUi and T.nlinj = djUjF, where T.F is the 
scalar field. 

• void Compute jrlin_RB(IncSF& T): same as Com- 
pute_nlin_scalar(IncSF& T). 

• void Compute jrhn(IncVF& W): nlin_i 

= (djUjUi — djWjWi) and W.nlin_i = 

(djUjWi — djWjUi), where Wi is the vector field 
associated with the Inc VF class W. 

• void Compute_nlin(IncVF& W, IncSFfe T): 

nlin_i = {djUjUi — djWjWi), and W.nlinJ = 

(djUjWi — djWjUi), and T.nlin = djUjF, with the 
same interpretation as given above. 

Class IncFluid: Incompressible Fluid 

The class IncFluid inherits IncVF and Time. Major 
functions of this class deal with time advancement of 
solver, forcing function, and input/output. The files and 
their associated functions are defined within this class. 
At present, the code includes Euler, Runge-Kutta sec- 
ond order (RK2), and Runge-Kutta fourth order (RK4) 
for the time advance function. The forcing function is 
used to include the buoyancy term in Rayleigh-Benard 
convection (RBC), Coriolis force in the rotating turbu- 
lence etc. 

For input /output, we have the option of reading/ writ- 
ing the data either the ASCII format or in High Density 
Format(HDF5) format. The HDF5 part of the code is 
being integrated with the main code. Also note that the 
classes IncVF and IncFluid have multiple inheritance. 

Solvers of TARANG 

We invoke the library functions discussed above to 
create solvers for fluid, magnetohydrodynamics, passive 
scalar, RBC flows etc. We illustrate a code segment con- 
taining the time-loop of fluid solver for an illustration. 

// A code segment of the fluid solver 
// Read initial condition 

U . Read_init_cond () ; 

int iter=0; // iterations 

U.Tnow = U.Tinit; 

do 

{ 

U . Compute _f orce () ; 
U . Compute_nlin () ; 
U . Add_f orce () ; 
U . Compute. pressure () ; 
U . Tdt = U . Get_dt () ; 
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U.Tnow = U.Tnow + U . Tdt ; 
iter ++ ; 

U. Time_ advance () ; 
// FIELD AT new time 
U . Output_all_inloop() ; 

} 

while (U.Tnow < U.Tfinal); 

Listing 2. Fluid Solver 

In the above code segment, U is an instantiation of the 
class IncFluid which contains the incompressible velocity 
field. Most of the functions are obvious. The function 
U.Get-dt() computes dt using CFL condition. 

For the RBC, the above code segment is modified 
slightly. We create an instantiation T of the class IncSF 
to represent the temperature field. 



// A code segment of the RBC solver 
// Read initial condition 

U . Read_init_cond (T) ; 

int iter=0; // iterations 

U.Tnow = U.Tinit; 

do 

{ 

U . Compute _f or ce (T) ; 

U . Compute_nlin (T) ; 

U . Add_f orce (T) ; 

U . Comput e_pressure () ; 

U . Tdt = U . Get_dt (T) ; 

U.Tnow = U.Tnow + U . Tdt ; 

iter ++ ; 

U. Time_advance (T) ; 
// FIELD AT new time 
U. Output_all_inloop (T) ; 

} 

while (U.Tnow < U.Tfinal); 



Listing 3. RBC solver 



SAMPLE SIMULATION RESULTS AND 
VALIDATION 

We have performed simulations on fluids, convective, 
and MHD flows using TARANG on grids from 64 3 to 
1024 3 [3H7]. The reader is referred to the published work 
for the scientific details. In the following discussion we 
detail some of the validations we performed before we 
launched large simulations. 



Simulations of RBC 

Rayleigh-Benard convection (RBC) is an idealized ver- 
sion of the thermal convection in fluid. In the set up, 
a layer of incompressible fluid is confined between two 



thermally conducting plates separated by a distance d. 
The bottom plate is heated and an adverse temperature 
gradient f3 is set across the fluid layer. The system is 
governed by the following equations: 

d t u + (u • V)u = -Vp + PRBz + PV 2 u, (1) 
d t 9 + (u • V)0 = u 3 + V 2 6>, (2) 
V ■ u = 0, (3) 

where u = (ui, 1^,113) is the velocity field, 9 is the per- 
turbation in the temperature field from the steady con- 
duction profile, and z is the vertically directed unit vec- 
tor. The equations are nondimensionalized by choosing 
length scale as d, velocity scale as d 2 /k, and temperature 
scale as j3d, where k is the thermal diffusivity of the fluid. 
Two non dimensional parameters in the equations are the 
Rayleigh number, R — ag/3d /un and the Prandtl num- 
ber, P — v/k, where a is the coefficient of the volume 
expansion, g is the acceleration due to gravity, and v is 
the kinematic viscosity of the fluid . We also use another 
parameter, reduced Rayleigh number r = R/R c , where 
R c is the critical Rayleigh number. 

The top and bottom boundaries are considered to be 
stress free and perfectly conducting: 

«3 = d z ui — d z u 2 =8 = 0, at z = 0,1. (4) 

We assume periodic boundary conditions along the hori- 
zontal direction (x and y-direction) . The boundary con- 
ditions chosen here are ideal. However they allow us to 
choose Fourier and sin or cos basis functions (SCFT) in 
our DNS. 

The amount of heat transported in the convection pro- 
cess is measured by the Nusselt number (Nu), which is 
defined as the ratio of total heat flux to the conductive 
heat flux. Using the nondimensionalization defined ear- 
lier, it can be shown 



Nu = l + (u 3 e), 



(5) 



where, (•) stands for spatial averaging. Please refer to 
Thual [5] , for a detailed derivation of the expression for 
Nusselt number. 

Thual [5] numerically solved the Eqs. l|3 under the 
above boundary conditions (Eq. |4| using a pseudo- 
spectral code. His simulations were performed in a two- 
dimensional (2D) box (aspect ratio T = 2^/2) for r = 1.1 
to 70 and P = 6.8. To verify TARANG we compare 
Nusselt numbers obtained in our simulations with those 
of Thual's. The above set of Eqs. (|T][3]) are solved numer- 
ically using TARANG under the above boundary condi- 
tions (Eq. [4]). We use same geometry (i.e. 2D box with 
aspect ratio 2y2) as used by Thual. We use Fourier basis 
functions for representation along the x, and sin or cos 
functions for representation along the z direction (SCFT 
basis). The validation results are shown in Table [I] Nu 
values obtained with TARANG are in very good agree- 
ment with those of Thual's until oscillation sets in the 
system. 
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TABLE I. Validation of TARANG against Thual's |8] 2D sim- 
ulations. We compare Nusselt numbers (Nu) obtained in 
our simulations with grid resolution 64 2 against Thual's 16 2 
(THU1), 32 2 (THU2), and 64 2 (THU3) simulations. All Nu 
values tabulated below are for P = 6.8. 

r THU1 THU2 THU3 TARANG 
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2.142 






2.142 


3 


2.678 






2.678 


4 


3.040 


3.040 




3.040 
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3.553 


3.553 




3.553 


10 


4.247 


4.244 




4.243 


20 


5.363 


5.333 


5.333 


5.333 


30 


6.173 


6.105 


6.105 


6.105 


40 


6.848 


6.742 


6.740 


6.740 


50 


7.441 


7.298 


7.295 


7.295 


70 




oscil. 


oscil. 


8.267 



PARALLELIZATION 

TARANG has been organized in a modular manner, 
so parallelization of the code was quite straight for- 
ward. Another major advantage was availability of par- 
allel FFTW. We essentially adopt FFTW's strategy for 
dividing the arrays etc. If p is the number of available 
processors, we divide each of the arrays into p segments. 
For example, a complex array A(N\,N2, Ns/2+1) is split 
into A(Ni/p, N 2 , N3/2 + 1) segments, each of which is 
handled by a processor. The other major parallel tasks 
needed is the multiplication in real space, which is han- 
dled by individual processors. Input/output is presently 
handled by the master node that collects/distributes data 
from the processor nodes. We are planning to implement 
parallel input /output using HDF5 functions. 

PORTING TO GPU 

Graphics Processing Units (GPUs) are getting popular 
in high performance computing due to their larger num- 
ber of cores. We have ported the FFT part of TARANG 
to GPUs, and have observed a reasonable speedup. We 
are in the process of porting the entire code to multiple 
GPU platform. 

FUTURE PLANS 

We have successfully performed simulations for fluid, 
MHD, and convective flows for grids up to 1024 3 on 



several platforms including PARAM YUVA (Centre for 
Advanced Computing, Pune), EKA (Computational Re- 
search Laboratory, Pune), HPC and CHAOS (both at 
IITK Kanpur). Attempts are being made to run turbu- 
lence simulations on higher grids. 

We have solved channel flow using Chebyshev and 
Fourier basis functions. We will be porting the full imple- 
mentation of the above basis function to be able to solve 
RBC and MHD flows under no-slip or mixed boundary 
conditions at the walls. The other planned modules are 
flows for the cylindrical and spherical geometry. We are 
also attempting to test and operationalize the magneto- 
convection module. 



CONCLUSIONS 



TARANG exploits the object-oriented programming 
features of C++ to build flow solvers for incompressible 
fluid, MHD, and convection. We adopt a modular ap- 
proach where general purpose functions are assembled to 
create solvers for different situations. This approach is 
proving to be very useful for constructing a large scale 
parallel softwares for fluid flows. 
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